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ABSTRACT 

We present high-resolution mass reconstructions for five massive cluster-lenses 
spanning a redshift range from z — 0.18-0.57 utilising archival Hubble Space Telescope 
(HST) data and applying galaxy-galaxy lensing techniques. These detailed mass mod- 
els were obtained from the observations by combining constraints from the strong and 
weak lensing regimes. We ascribe local weak distortions in the shear maps to pertur- 
bations induced by the presence of galaxy haloes around individual bright early-type 
cluster member galaxies. This technique constrains the mass enclosed within an aper- 
ture for these subhaloes. We are sensitive to a specific mass range for these subhaloes, 
10 11 - 1O 12 ' 5 M , which we associate with galaxy-scale subhaloes. Adopting a para- 
metric model for the subhaloes, we also derive their velocity dispersion function and 
the aperture radius function. The mass spectrum of substructure in the inner regions 
of the observed clusters is directly compared with that in simulated clusters extracted 
from the Millennium Simulation. The mass function, aperture radii and velocity dis- 
persion function are compared in detail. Overall, we find good agreement between the 
distribution of substructure properties retrieved using the lensing analysis and those 
obtained from the simulation. We find that the fraction of total cluster mass asso- 
ciated with individual subhaloes within the inner 0.5 - 0.8ft. _1 Mpc of our clusters 
ranges from 10-20 per cent, in broad agreement with simulations. Our work provides 
a powerful test of the ACDM model, which appears consistent with the amount of 
observed substructure in massive, lensing clusters based on present data. 

Key words: gravitational lensing, galaxies: fundamental parameters, haloes, meth- 
ods: numerical 



1 INTRODUCTION 

Gravitational lensing has emerged as one of the most power- 
ful techniques to map mass distributions on a range of scales: 
galaxies, clusters and beyond. The distortion in the shapes 
of background galaxies viewed through foreground mass dis- 
tributions is independent of the dynamical state of the lens, 
therefore, compared with other methods for mass estima- 
tion there are fewer biases in lensing mass determinations. 
Here, we focus on mapping in detail the mass distribution 
inside the inner regions of massive clusters of galaxies us- 
ing Hubble Space Telescope (HST) observations. We exploit 
the technique of galaxy-galaxy lensing, which was originally 
proposed as a method to constrain the masses and spatial 
extents of field galaxies (Brainerd, Blandford & Smail 1996), 
which we have since extended and developed over the years 
to apply inside clusters (Natarajan & Kneib 1997; Natarajan 
et al. 1998; 2002a). 

The detailed mass distribution within clusters and 



specifically the fraction of the total cluster mass associated 
with individual galaxies has important implications for the 
frequency and nature of galaxy interactions in clusters (Mer- 
ritt 1983; Richstone 1976; Farouki & Shapiro 1981; Moore et 
al. 1996; Ghigna et al. 1998; Okamato & Habe 1999). Knowl- 
edge of the dynamical history of clusters enables a deeper 
understanding of the physical processes that shape their as- 
sembly and evolution. The discovery of strong evolution be- 
tween z ~ 0.5 and the present-day in the morphological (and 
star-formation) properties of the galaxy populations in clus- 
ters has focused interest on environmental processes which 
could effect the gaseous component and dark matter halo of 
a cluster galaxy (e.g. Couch et al. 1994, 1998). 

The global tidal field of a massive, dense cluster poten- 
tial well is expected to be strong enough to truncate the dark 
matter halo of a galaxy whose orbit penetrates the cluster 
core. Therefore, probing the extents of galaxy haloes in clus- 
ters can provide invaluable clues to dynamically dominant 
processes in clusters. For instance, the survival of individual, 
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compact dark haloes associated with cluster galaxies sug- 
gests a high probability for galaxy-galaxy collisions within 
rich clusters over a Hubble time. However, since the inter- 
nal velocity dispersions of cluster galaxies ( £ 200 km s~ 1 ) are 
significantly lower than their orbital velocities, these interac- 
tions are, in general, unlikely to lead to mergers, but rather 
encounters of the kind simulated in the galaxy harassment 
picture by Moore et al. (1996, 1998). 

Previous work on galaxy-galaxy lensing in the moderate 
redshift field has identified a signal associated with massive 
haloes around typical field galaxies, extending to beyonqj 
100 kpc (e.g. Brainerd, Blandford & Smail 1996; Ebbels et 
al. 2000; Hudson et al. 1998; Hoekstra et al. 2004). In par- 
ticular, Hoekstra et al. (2004) report the detection of a finite 
truncation radius of 185 ± 30 kpc via weak lensing by galax- 
ies based on imaging data from the Red-Sequence Cluster 
Survey. Galaxy-galaxy lensing results from the analysis of 
the Sloan Digital Sky Survey data (McKay et al. 2002; Shel- 
don et al. 2004; Guzik & Seljak 2002) have contributed to a 
deeper understanding of the relation between mass and light. 
Similar analysis of galaxies in the cores of rich clusters sug- 
gests that the average mass-to-light ratio and spatial extent 
of the dark matter haloes associated with morphologically 
classified early-type galaxies in these regions may differ from 
those of field galaxies with comparable luminosity (Natara- 
jan et al. 1998, 2002a). We find that at a given luminosity, 
galaxies in clusters have more compact halo sizes and lower 
masses (by a factor of 2-5) compared to their field counter- 
parts. The mass-to-light ratios inferred for cluster galaxies 
in the V-band are also lower than those of field galaxies 
with comparable luminosity. This is a strong indication of 
the effect of the dense environment on the properties of dark 
matter haloes. 

In this paper, we present a determination of the mass 
function and other detailed properties of substructure in 
clusters using galaxy-galaxy lensing techniques. A high res- 
olution mass model tightly constrained by strong and weak 
lensing observations is constructed including individual clus- 
ter galaxies and their associated dark matter haloes. These 
lensing models are new and incorporate recent observed 
spectroscopic redshifts of several additional multiple images. 
We show that over a limited mass range we can successfully 
construct the mass function of subhaloes inside clusters. In 
earlier work, we compared the mass function obtained from 
lensing with that extracted from a massive cluster simulated 
at extremely high resolution (Natarajan & Springel 2004). 
In this work, we are able to compare with a large ensemble of 
simulated clusters, significantly strengthening the statistical 
significance of our results. 

N-body simulations represent an indispensable tool 
for investigating the non-linear growth of structures in its 
full geometrical complexity. The high numerical resolution 
achieved in recent years has demonstrated that the cores 
of dark matter haloes that fall into a larger system can 
survive for a relatively long time as self-gravitating objects 
orbiting in the smooth dark matter background of the sys- 
tem. A wealth of dark matter substructures is now routinely 
detected and studied with the aid of N-body simulations 



1 We adopt h = flo/lOOkms -1 Mpc -1 = 0.7 and H A = 0.7, and 
scale other published results to this choice of parameters. 



confirming that the existence of substructures is a generic 
prediction of hierarchical structure formation in Cold Dark 
Matter (CDM) models. 

The subhalo mass function (i.e. the abundance of dark 
matter substructure as a function of mass) represents an 
important prediction of hierarchical CDM structure forma- 
tion models and has been subject of intense studies since 
the 'dwarf galaxy crisis' was identified (Moore et al. 1999; 
Klypin et al. 1999; Stoehr et al. 2003). Within a radius 
of 400 ft -1 kpc, from the Milky Way, cosmological models 
of structure formation predict ~ 30 dark matter satellites 
with circular velocities in excess of 20kms _1 and mass 
greater than 3 x 10 8 Mq. This number is significantly higher 
than the dozen or so satellites actually observed around 
our Galaxy. Different explanations have been suggested for 
this discrepancy. The missing satellites could for example be 
identified with the detected High Velocity Clouds (Mailer & 
Bullock 2004). Warm or self -interacting dark matter could 
also selectively suppress power on the small scales, therefore 
reducing the predicted number of satellites (Spergel & Stein- 
hardt 2000). The leading hypothesis however remains that 
the solution to this problem lies in astrophysical processes 
such as heating by a photo-ionising background that sup- 
presses star formation in small haloes at early times (Bullock 
et al. 2000; Benson et al. 2002; Kravtsov et a. 2004). On the 
scale of galaxy clusters, many more dark matter structures 
are expected to host visible galaxies, thus making the com- 
parison with expectations from numerical simulations less 
affected by uncertainties in the physics of galaxy formation. 
Full consistency is to be expected here, therefore providing 
an important test of the CDM paradigm. 

Our innovative application of gravitational lensing en- 
ables a mapping of the distribution in mass of dark matter 
substructure (subhalo mass function) therefore allowing a 
direct comparison with results from numerical simulations. 
The strength of the lensing analysis presented here derives 
from the combination of both strong and weak lensing fea- 
tures which are used together to construct a high resolution 
mass map of a galaxy cluster. Anisotropics in the shear field 
(i.e. departure from the coherent tangential signal) in the 
vicinity of bright, early- type cluster members are attributed 
to the presence of these local potential wells. Statistically 
stacking this signal provides a way to quantify the masses as- 
sociated with individual galaxy haloes. This is accomplished 
using a maximum likelihood estimator to retrieve character- 
istic properties for a typical subhalo in the cluster. 

The comparison with clusters from the Millennium Sim- 
ulation we use here is particular powerful due to the large 
volume and high resolution of this simulation. The simula- 
tion used more than 10 billion particles to trace the evolution 
of the dark matter distribution in a cubic region of the Uni- 
verse over 2 billion light-years on a side. This makes it an 
ideal data set for comparison with lensing clusters, provid- 
ing dozens of highly resolved, massive lensing clusters out 
to redshifts z ~ 0.6 and beyond. 

The outline of this paper is as follows: in Section 2, we 
describe briefly the formalism for analysing galaxy-galaxy 
lensing in observed clusters including a synopsis of the 
adopted models. Section 3 presents the best-fit lens mod- 
els and discusses the uncertainties and sources of error. A 
detailed comparison with clusters from the Millennium Sim- 
ulation is presented in Section 4. The results of comparing 
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the mass function of substructure, the distributions of aper- 
ture radii and the velocity dispersions are discussed. Finally, 
we conclude with a discussion of the implications of our re- 
sults for the ACDM model and the future prospects of this 
work in Section 5. 



2 GALAXY-GALAXY LENSING IN CLUSTERS 

2.1 Framework for analysis 

In this section, we briefly outline the analysis framework, 
noting that further details can be found in earlier papers 
(Natarajan & Kneib 1997; Natarajan et al. 1998). For the 
purpose of extracting the properties of the subhalo popula- 
tion in clusters, a range of mass scales is modelled parametri- 
cally. The X-ray surface brightness maps of clusters suggest 
the presence of a smooth, dominant, large scale mass com- 
ponent. Clusters are therefore modelled as a super-position 
of a smooth large-scale potential and smaller scale potentials 
that are associated with bright early-type cluster members: 



0smooth ~T ^i ' 



where <^ sm0 oth is the potential of the smooth component and 
4> Pi is the potential of the subhalo associated with the ith 
galaxy, and is treated as perturber. The amplification matrix 
A^ 1 can be decomposed into contributions from the main 
clump and the perturbcrs: 



A' 1 = (1 



^smooth 2-<i^p,M 7: 



smooth -J20 



hJlt 



where k is the magnification and 7 the shear. The shear 7 is 
in fact a complex number and is used to define the reduced 
shear <;, which is the quantity that is measured directly from 
the observed shapes of background galaxies. Similarly, the 
reduced shear can also be decomposed as 



Stot 



7 



/smooth ' % 'Pi 



In the frame of an individual perturber j (and neglecting 
the effect of perturber i if i ^ j), this simplifies to: 



/smooth 1 Jpj 



J- ^smooth K; 



pj 



Restricting our analysis to the weak regime, and thereby 
retaining only the first order terms from the lensing equation 
for the shape parameters (e.g. Kneib et al. 1996), we have: 



7 c e l c for the smooth cluster component over the area of in- 
tegration. In the frame of the perturber, the averaging pro- 
cedure allows efficient subtraction of the large-scale compo- 
nent, enabling the extraction of the shear component in- 
duced in the background galaxies only by the local per- 
turber. The background galaxies are assumed to have in- 
trinsic ellipticities drawn from a known distribution (see the 
next section for further details). Schematically, the effect of 
the cluster on the intrinsic ellipticity distribution of back- 
ground sources is to cause a coherent displacement r and 
the presence of perturbers merely adds small-scale noise to 
the observed ellipticity distribution. 

The feasibility and robustness of signal detection has 
been amply demonstrated in earlier papers by Natarajan 
et al. The primary limitations in this analysis arise from 
the total number of distorted background galaxies, and the 
accuracy with which the smooth cluster component can 
be constrained using multiple images in the inner region. 
The partition into this smooth component and its effective 
subtraction in fact boosts the shear induced by the per- 
turber. In particular, the shear induced by the subhalo has 
a (^smooth + K P A term in the denominator, which becomes 
non-negligible in the cluster centre. The subtraction of the 
larger-scale component reduces the noise in the polarisa- 
tion measure by about a factor of two in cluster cores. This 
differenced averaging prescription for extracting the distor- 
tions induced by the possible presence of dark haloes around 
cluster galaxies is feasible with HST quality data as we have 
shown in earlier work (Natarajan et al. 1998, 2002a). The 
^i 7pi ^29 p . , robustness of this technique also has been demonstrated in 
our earlier published work. Note here that it is the pres- 
ence of the underlying large-scale smooth mass distribution 
(with a high value of « c ) that enables the extraction of the 
weak signal riding on it. It is instructive to keep in mind 
that, in the regimes of interest discussed here, the distor- 
tion induced by the cluster-scale smooth component for a 
pseudo-isothermal elliptical component model (PIEMD) in 
the inner-most (with a velocity dispersion of 1000 kms -1 and 
at R/rt < 0.1) regions is typically of the order of 20 - 40 per 
cent or so in background galaxy shapes, and the perturbers 
produce distortions (smaller scale PIEMDs with a velocity 
dispersion of 220 kms -1 ) of the order of 5 - 10 per cent, 
significantly more than in the case of weak-lensing by large 
scale structure or cosmic shear, wherein the distortions are 
of the order of 1 per cent. 



(1) 



(2) 



(3) 



9i = 9s +9tot, ( 4 ) 

where TjJ 2 ] is the distortion of the image, ~g s the intrinsic 
shape of the source, and g tot is the distortion induced by 
the lensing potentials. 

In the local frame of reference of the perturbers, the 
mean value of the quantity g 7 and its dispersion can be com- 
puted in circular annuli (at radius r from the perturber cen- 
tre) strictly in the weak-regime, assuming a constant value 



2 The measured image shape and orientation are used to con- 
struct a complex number whose magnitude is given in terms of 
the semi-major axis (a) and semi-minor axis (b) of the image and 
the orientation is the phase of the complex number. 



2.2 Modelling the cluster 

Each of the clusters studied in this paper preferentially 
probes the high mass end of the cluster mass function and 
has a surface mass density in the inner regions which is 
higher than the critical value, therefore producing a number 
of multiple images of background sources. By definition, the 
critical surface mass density for strong lensing is given by: 



^crit 



D s 



(•») 



4ttG D d D ds ' 

where D 3 is the angular diameter distance between the ob- 
server and the source, Dd the angular diameter distance be- 
tween the observer and the deflecting lens, and Dd s the 
angular diameter distance between the deflector and the 
source. When the surface mass density in the cluster is in 
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excess of this critical value, strong lensing phenomena with 
high magnification are observed. 

In general two types of lensing effects are produced - 
strong: multiple images and highly distorted arcs; and weak: 
small distortions in background image shapes determined 
by the criticality of the region. Viewed through the central, 
dense core region of the mass distribution, where k > 1 
strongly lensed features are observed. Note that the inte- 
grated lensing signal detected is due to all the mass dis- 
tributed along the line of sight in a cylinder projected onto 
the lens plane. In this and all other cluster lensing work, the 
assumption is made that individual clusters dominate the 
lensing signal as the probability of encountering two mas- 
sive rich clusters along the same line-of-sight is extremely 
small due to the fact that these are very rare objects in 
hierarchical structure formation models. 

With our current sensitivity limits, galaxy-galaxy lens- 
ing within the cluster is primarily a tool to determine the 
total enclosed mass within an aperture. We lack sufficient 
sensitivity to constrain the detailed mass profile for indi- 
vidual cluster galaxies. With higher resolution data in the 
near future we will be able to obtain constraints on the 
slopes of mass profiles in subhaloes. In this paper, we there- 
fore concentrate on pseudo-isothermal elliptical components 
(PIEMD models, derived by Kassiola & Kovner 1993) appro- 
priately scaled for both the main cluster and the substruc- 
tures. Some of the cluster studied in this work are bi-modal 
(i.e. have two significant, large-scale mass peaks), therefore, 
when required we employ two large-scale smooth potentials 
and a super-position of subhaloes. We find that the results 
obtained for the characteristics of the subhaloes (or per- 
turbers) is largely independent of the form of the mass dis- 
tribution (the aperture mass is comparable for the NFW 
and PIEMD models) used to model the smooth, large-scale 
component. 

To quantify the lensing distortion induced by the global 
potential, both the smooth and individual galaxy-scale 
haloes are modelled self-similarly using a surface density 
profile, £(.R), which is a linear superposition of two PIEMD 
distributions, 



E(R) = 



S r 



r o/rt I ^r 2 + R 2 ^Jr 2 + R 2 



(6) 



with a model core-radius ro and a truncation radius r t ^> 
ro- These parameters (r t ,ro) are tuned for both the smooth 
component and the perturbers to obtain mass distributions 
on the relevant scales. The coordinate R is a function of x, 
y and the ellipticity, 



R* 



+ 



y 



(l + e) 2 (1-6) 



where 



a + b 



(7) 



The mass enclosed within radius R for the e — model is 
given by 



M(R) 



2nT,oro 



\A§ + R 2 - y/ r t + R 2 + (n - ro) 



.(8) 



One of the attractive features of this model is that the total 
mass M is finite, M oc Soron. Besides, analytic expressions 
can be obtained for the all the quantities of interest, k, 7 



and g, e.g. 



k(R) — Kq 



TO 



(1- ro/n) [y/rf+R 2 y/rf+R 2 



2/to 



4ttG Di s D, 



Is-LAdI 



D 



(9) 



(10) 



where D\ s , D os and D \ are the lens-source, observer- source 
and observer-lens angular diameter distances, respectively, 
which do depend on the choice of cosmological parameters. 
To obtain g(R) knowing the magnification k(R), we solve 
Laplace's equation for the projected potential 02D, evalu- 
ate the components of the amplification matrix and then 
proceed to solve directly for y(R), and then g(R): 



■y(R) = k 



+ 



/R^TTf R 2 



+ 4(\A R2 + r o-n>) 



VR T +^ 



R 2 



( V / R 2 +r 2 -r t ) 



Scaling this relation by rt gives for ro < R < rt: 



l(R/rt) oc 



So r t 



a 

7F 



t)-lR 

where a is the velocity dispersion and for r ( ) < r t < R 



■y(R/r t ) oc 



Ep_r£; 
77 R? 



Mu 



/?- 



(11) 



(12) 



(13) 



where M to t is the total mass. In the limit that R ^> rt, we 
have 



7 (fl) 



3kq . 2 
2~R3 [r ° 



2, . 2k 



00, f(R) 



- ro], (14) 

0, g(R) -* and t(R) -> 0, as 



and as R 

expected. 

It is further assumed that the ellipticity and the ori- 
entation of the dark matter subhaloes associated with the 
early-type cluster members is identical to that of the galax- 
ies themselves. Additionally, in order to relate the light dis- 
tribution to key parameters of the mass model above, we 
adopt a set of physically motivated scaling laws for the clus- 
ter galaxies (Brainerd et al. 1996): 



(To = (T()» 



rtAjz) -(I 5 ) 



These in turn imply the following scaling for the r t /ro ratio 
r,: 

r r , \L*J 

The total mass M ap enclosed within an aperture r t » and the 
total mass-to- light ratio M/L then scale with the luminosity 
as follows: 



Wa P oc (ToVt, (77) , M/L oc (To«r t , ( — 



L* 



a-1/2 



(17) 



where a determines the size of the galaxy halo. For a — 
0.5 the assumed galaxy model has constant M/L with lu- 
minosity (but not as a function of radius) for each galaxy. 
We adopt a = 0.5 throughout this work, as this scaling 
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is empirically motivated by the Faber- Jackson relation for 
early- type galaxies (Brainerd, Blandford & Smail 1996). We 
assume these scaling relations and recognise that this could 
ultimately be a limitation of our model but the evidence 
at hand supports the fact that mass traces light efficiently 
both on cluster scales (Kneib et al. 2003) and on galaxy 
scales (McKay et al. 2002; Wilson et al. 2001). The calibra- 
tion factors Co* and r t , in eqn. (15) are provided by the 
results of a likelihood method discussed further in Section 
2.3. 



2.2.1 The intrinsic shape distribution of background 
galaxies 

As in all lensing work, it is assumed here as well that the 
intrinsic or undistorted distribution of shapes of background 
galaxies is known. This distribution is obtained from shape 
measurements taken from deep images of blank field surveys. 
Previous analysis of deep survey data such as the MDS fields 
(Griffiths et al. 1994) showed that the ellipticity distribu- 
tion of sources is a strong function of the sizes of individual 
galaxies as well as of their magnitude (Kneib et al. 1996). 
For the purposes of our modelling, the intrinsic ellipticities 
for background galaxies are assigned in concordance with an 
ellipticity distribution p(t$) where the shape parameter r is 
defined as r = (a 2 — b 2 )/(2ab) derived from the observed 
ellipticities of the CFHT12k data (see Limousin et al. 2005 
for details): 



For the normalised redshift distribution at a given mag- 
nitude m (in the given band) we therefore have 



p(rs) = i~s exp 



T.s- 



v = 1.15, 8 = 0.25. 



(18) 



Note that this distribution includes accurately measured 
shapes of galaxies of all morphological types. In the like- 
lihood analysis this distribution p(rs) is the assumed prior, 
which is used to compare with the observed shapes once the 
effects of the assumed mass model are removed from the 
background images. We note here that the exact shape of 
the ellipticity distribution, i.e. the functional form and the 
value of 8 and v do not change the results, but alter the con- 
fidence levels we obtain. The width of the intrinsic ellipticity 
distribution on the other hand is the fundamental limiting 
factor in the accuracy of all lensing measurements. 



2.2.2 The redshift distribution of background galaxies 

While the shapes of lensed background galaxies can be mea- 
sured directly and reliably by extracting the second moment 
of the light distribution, the precise redshift for each weakly 
object is in general unknown and therefore needs to be as- 
sumed. Using multi-waveband data from surveys such as 
COMBO-17 (Wolf et al. 2004), photometric redshift esti- 
mates can be obtained for every background object. Typi- 
cally the redshift distribution of background galaxies is mod- 
elled as a function of observed magnitude P(z, m). We have 
used data from the high-redshift survey VIMOS VLT Deep 
Survey (Le Fevre et al. 2004) as well as recent CFHT12k R- 
band data to define the number counts of galaxies, and the 
HDF prescription for the mean redshift per magnitude bin, 
and find that the simple parameterisation of the redshift dis- 
tribution used by Brainerd, Blandford & Smail (1996) still 
provides a good description to the data. 



N{z)\m = 
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/?(* )exp[-(^ 
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= 1.5 and 
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L dmjj 


- m R0 ) 



(19) 



(20) 



with Zmedian being the median redshift, and dz mc dian/dm_R 
being the change in median redshift with say the _R-band 
magnitude, ma- 

However, we note here in agreement with another recent 
study of galaxy-galaxy lensing in the field by Kleinheinrich 
et al. (2005), that the final results for the aperture mass 
presented here are primarily sensitive to the choice of the 
median redshift of the distribution rather than the individ- 
ual assigned values. 



2.3 The maximum-likelihood method 

There are two key aspects to constructing a successful lens 
model for the clusters analyzed here. Firstly, we must iden- 
tify multiple-imaged background sources with reliable red- 
shift measurements whose properties can be used to con- 
strain the total projected mass within the Einstein radius. 
The strong lensing observations provide constraints in the 
inner regions for the global smooth potential. Secondly, we 
use the weak shear detected out to larger radii to constrain 
the behavior of the smooth component and simultaneously 
statistically constrain the properties of the subhaloes. This 
is done by quantifying the anisotropies that arise in the shear 
field and attributing these anisotropies to the presence finite 
mass subhaloes. 

Parameters that characterise both the global compo- 
nent and the perturbers are optimised, using the observed 
strong lensing features - positions, magnitudes, geometry of 
multiple images and measured spectroscopic redshifts, when 
known, along with the smoothed shear field as constraints. 
As initial values of the mass model we provide a center, 
velocity dispersion, ellipticity, orientation, and truncation 
radius for the main clump with tolerance ranges. For the 
subhaloes, we provide the locations (coincident with the po- 
sitions of the selected early-type cluster members) , elliptic- 
ity and orientation (taken from measurements for the cluster 
early- types). In combination with the parameterisation pre- 
sented in the previous section, we then optimise and extract 
values for the central velocity dispersion and the aperture 
scale (ao»,r t ») for a typical L* cluster galaxy. 

A maximum-likelihood estimator is used to obtain sig- 
nificance bounds on fiducial parameters that characterise 
a typical L* subhalo in the cluster. We have extended the 
prescription proposed by Schneider & Rix (1997) for galaxy- 
galaxy lensing in the field to the case of lensing by galaxy 
subhaloes in the cluster (Natarajan & Kneib 1997, Natara- 
jan et al 1998). The likelihood function of the estimated 
probability distribution of the source ellipticities is max- 
imised for a set of model parameters, given a functional 
form of the intrinsic ellipticity distribution measured for 
faint galaxies. For each 'faint' galaxy j, with measured shape 
T bs, the intrinsic shape rg. is estimated in the weak regime 
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by subtracting the lensing distortion induced by the smooth 
cluster model and the galaxy subhaloes, 



TSi — T bs 



■^ N c 



7c 



(21) 



where S^ ■y Pi is the sum of the shear contribution at a given 
position j from N c perturbers. This entire inversion proce- 
dure is performed numerically using code that builds on the 
ray-tracing routine lenstool written by Kneib (1993). This 
machinery accurately takes into account the non-linearities 
arising in the strong lensing regime. Using a well-determined 
'strong lensing' model for the inner regions of the clus- 
ters derived from the positions, shapes and magnitudes of 
the highly distorted multiple-imaged objects along with the 
shear field determined from the shapes of the weakly dis- 
torted background galaxies and assuming a known func- 
tional form for p(rs), the probability distribution for the 
intrinsic shape distribution of galaxies in the field, the like- 
lihood for a guessed model is given by 



C(a 0t ,r tf ) = n^p(r Sj ), 



(22) 



where the marginalisation is done over (cto*, ?"**)• We com- 
pute C assigning the median redshift corresponding to the 
observed source magnitude for each arclet. The best fitting 
model parameters are then obtained by maximising the log- 
likelihood function with respect to the parameters ao* and 
r t ». Note that the parameters that characterise the smooth 
component are also simultaneously optimised. The likeli- 
hood can also be marginalised over a complementary pair 
of parameters, e.g. using the luminosity scaling index a and 
the aperture mass M ap directly. In this work, we explore 
both choices. 



3 BEST-FIT LENSING MASS MODELS 

A composite mass model is constructed for the clusters start- 
ing with the super-posed PIEMDs. The strong lensing data, 
i.e. the geometry, positions, relative brightness, redshifts and 
parities of the multiple images are used to obtain the mass 
enclosed within the Einstein radius which is used as an initial 
constraint for the integrated mass in the inner regions. The 
contribution to the shear and magnification from all poten- 
tials (large-scale and small-scale perturbers) is calculated at 
the location of every observed background source galaxy and 
the inversion of the lensing equation is performed. The ob- 
served shape and magnification of each and every distorted 
background galaxy is compared to that computed from the 
model and the subhalo mass distribution is modified itera- 
tively until the best match between the observations and the 
model is found simultaneously for all background sources. 
Strong lensing constraints principally drive the likelihood 
results. Therefore, we use many sets of multiple images with 
measured redshifts for each cluster as inputs. Compared to 
our earlier work (Natarajan et al. 2002) the models pre- 
sented here have several improvements, we have incorpo- 
rated more multiple image families with redshifts that have 
since become available and we have also modified the al- 
gorithm to enable finer grid searching. In our analysis, the 
center of the cluster is picked to coincide with the location 
of the brightest cluster galaxy. 

The basic steps of our analysis involve lens inver- 
sion, modeling and optimization, which are done using the 



lenstool software utilities (Kneib 1993). These utilities are 
used to perform the ray tracing from the image plane to 
the source plane with a specified intervening lens. This is 
achieved by solving the lens equation iteratively, taking into 
account the observed strong lensing features, positions, ge- 
ometries and magnitudes of the multiple images. In some 
cases, we also include a constraint on the location of the crit- 
ical line to tighten the optimization. Additionally, we fix the 
core radius of an L* subhalo to be 0.1 kpc, as by construc- 
tion our analysis cannot constrain this quantity. In addition 
to the likelihood contours, the reduced % 2 f° r the- best-fit 
model is also robust. We describe some pertinent features of 
each cluster and their respective mass models below. 
A 2218 

Our best fit mass model for the cluster is bimodal, com- 
posed of two large scale clumps around the cD and the sec- 
ond brightest cluster galaxy. This model is an updated ver- 
sion of that constructed by Kneib et al. (1996). It includes 
40 additional small-scale clumps that we associate with lu- 
minous early-type galaxies in the cluster core. Only about 
10 per cent of the total cluster mass is in substructures, i.e. 
associated with galaxy scale haloes. The aperture mass, in- 
tegrated over the truncation radius r ap = 40 kpc, yields a 
characteristic mass of 1.4 x 10 12 Mq, with a total mass-to- 
light ratio in the V-band of ~ 5. 8 ± 1.5 and a central velocity 
dispersion of about 180 km s -1 . 
A 2390 

The cluster has an unusual feature - a strongly lensed 
almost 'straight arc' (Pello et al. 1991) approximately 38 
arcsec (~ 170 kpc) away from the central galaxy, in addi- 
tion to many other arcs and arclets that have been utilised in 
our modeling. We find a best-fit mass model with two large- 
scale components, that yield a projected mass within the 
radius defined by the brightest arc of ~ 1.8 ±0.2 x 10 14 M . 
Our best-fit composite lensing model for A 2390 incorpo- 
rates 40 perturbers associated with early-type cluster mem- 
bers whose characteristic parameters are optimized in the 
maximum-likelihood analysis. The integrated mass within 
the ~ 18 kpc tidal radius for a typical L* cluster galaxy 
is about 6.4 x 10 1 M© giving a total mass-to-light ratio in 
the V-band of about 4.2 ± 1.3. Again, 90 per cent of the 
total mass of the cluster is consistent with being smoothly 
distributed. 
CI 2244-02 

The best-fit lens model for this cluster has two compo- 
nents which both have fairly low velocity dispersions. This is 
the least massive lensing cluster in the sample studied here. 
The X-ray mass estimate from the ASCA data (Ota et al. 
1998) is in good agreement with our best-fit lensing mass 
model. This is despite the fact that the X-ray temperature 
of CI 2244— 02 is at least a factor of two higher than that 
expected from the average luminosity-temperature relation. 

The tidal truncation radius obtained for a typical L* 
cluster galaxy in CI 2244 is the largest in the sample studied 
here and is 55 ± 12 kpc. This is consistent with the fact that 
the central density in CI 2244 is the lowest. The total mass- 
to-light ratio in the V-band for a fiducial L* is 3.2 ± 1.2. 
Approximately 20 per cent of the total mass is in substruc- 
ture within the mass range 10 11 — 10 12 ' 5 M©. 
CI 0024+16 

Our best fit mass model for the inner regions takes 
into account the small scale dark haloes associated with the 
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z 


a* 


(kpc) 


M^p/Ly 

(M©/L ) 


(10 n M Q ) 


^clus 

(kms -1 ) 


Pclus(r = 0) 
(10 6 M Q kpc" 3 ) 


A 2218 




0.17 


180 ± 10 


40 ±12 


5.8 ±1.5 


~ 14 


1070 ± 70 


3.95 


A 2390 




0.23 


200 ± 15 


18 ±5 


A 4.2 ±1.3 


~ 6.4 


1100 ±80 


16.95 


AC 114 




0.31 


192 ± 35 


17 ±5 


6.2 ±1.4 


~ 4.9 


950 ± 50 


9.12 


CI 2244- 


-02 


0.33 


110 ±7 


55 ±12 


3.2±1.2 


~ 6.8 


600 ± 80 


3.52 


Cl 0024+16 


0.39 


125 ±7 


45 ±5 


2.5 ±1.2 


~ 6.3 


1000 ± 70 


3.63 


CI 0054- 


-27 


0.57 


230 ± 18 


20 ±7 


5.2 ±1.4 


~ 9.4 


1100 ± 100 


15.84 



Table 1. Parameters that define the mass models of the subhalocs for the lensing clusters. 



early- type members in the core, and requires a two com- 
ponent model for the sub-clusters. Integrating the best-fit 
mass model shown, we find that (i) about 10 per cent of the 
total cluster mass is in galaxy-scale haloes and (ii) the total 
mass estimate is in good agreement with that obtained by 
Kneib et al. (2003) where data from a much larger field of 
view were used. 

Even on the large scales probed by Kneib et al. (2003) 
it was found that mass and light traced each other rather 
well at large radii. A typical L* cluster galaxy was found to 
have a truncation radius of 45 ±5 kpc, and a central velocity 
dispersion of 125 ± 7kms _1 . 
CI 0054 27 

The lensing signal from CI 0054—27 is best fit by a single 
smooth dark matter component and subhaloes associated 
with bright, early- type members making it the only uni- 
modal cluster in the sample studied here. The mass enclosed 
within ~ 400 kpc is of the order of 1.8 ± 0.4 x 10 14 0. 

The characteristic central velocity dispersion of a typi- 
cal L* galaxy in this cluster is higher than in A 2218, A 2390 
or C10024±16, all of which are by contrast bimodal in the 
mass distribution. In this cluster, about 20 per cent of the 
total mass is in substructure. However, CI 0054— 27 is the 
most distant cluster studied here and is likely to be still 
evolving and assembling, accounting for the high mass frac- 
tion in substructure. 

Note here that the choice of a determines only the scal- 
ing of the outer radius of a fiducial subhalo with luminosity. 
With the data used in this paper it is not possible to distin- 
guish between various values of a - some values are clearly 
more physical than others. Therefore, this implies that we 
are sensitive to the integrated mass within an aperture that 
is determined primarily by the anisotropy in the shear field 
and not by the details of how the subhalo masses are trun- 
cated. We also find that out to 500 kpc in all clusters only 
10-20 per cent of the total mass is associated with galaxy 
haloes. At this radius (to which we are limited due to the 
size of the HST WFPC2 fields) , most of the mass is in the 
large-scale component. This fraction is likely to be a strong 
function of cluster-centric radius. The dependence of these 
aperture radii with distance from the cluster can be explored 
with wide-field HST data and we are in the process of doing 
so for the cluster C10024±16 (Natarajan et al. 2006). 

3.1 Results from the lens models 

We successfully construct high resolution mass models for 
all five clusters from the unambiguous galaxy-galaxy lens- 



ing signal detected using the maximum-likelihood analysis. 
We constrain the mass enclosed within an aperture for a 
fiducial halo. The maximum-likelihood analysis yields the 
following: (i) the mass-to-light ratio in the V-band of a typ- 
ical L* galaxy does not evolve significantly as a function of 
redshift, (ii) the fiducial truncation radius of an L* galaxy 
varies from about 20 kpc to 70 kpc depending on the clus- 
ter (iii) the typical central velocity dispersion is roughly 
180 kms -1 . In Fig. 1, wc show the recovered parameters 
(cro*, rt*) for an L* galaxy from the maximum-likelihood 
analysis. For the galaxy mass model adopted in our analysis, 
the total mass of an L* galaxy varies from ~ 4.9 x 10 11 Mq 
to ~ 1.4 x 10 12 M . We note that the distribution of these 
recovered parameters points to roughly two kinds of typical 
haloes, compact, dense ones and more extended ones. How- 
ever, given that these lensing clusters do not constitute a 
well denned sample it is premature to claim any dichotomy. 
The mass-to-light ratios obtained take passive evolution of 
elliptical galaxies into account as given by the stellar popu- 
lation synthesis models of Bruzual & Chariot (2003), there- 
fore any detected trend reflects pure mass evolution. The 
mass obtained for a typical bright cluster galaxy by Tyson 
et al. (1998) using only strong lensing constraints inside the 
Einstein radius of the cluster C10024±1654, at z = 0.41, is 
consistent with our results. All error bars quoted here are 
~3cr. 

By construction, the maximum-likelihood technique 
presented here provides the mass spectrum of subhaloes in 
the cluster directly. Note that as stated before, in perform- 
ing the likelihood analysis to obtain characteristic parame- 
ters for the subhaloes in the cluster, it is assumed that light 
traces mass. This is an assumption that is well supported by 
galaxy-galaxy lensing studies in the field (Wilson et al. 2001) 
as well as in clusters (Clowe & Schneider 2002). In fact, all 
lens modeling and rotation curve measurements suggest an 
excess of baryons in the inner regions. Note however that 
for our choice of mass model (the PIEMD) the mass to light 
ratio is not constant with radius within an individual galaxy 
halo. Since the procedure involves a scaled, self-similar mass 
model that is parametric, we obtain a mass estimate for the 
dark haloes (subhaloes) as a function of the luminosity of 
the early-type galaxy hosted by them. This provides us with 
a clump mass spectrum. We surmise that truncation by the 
cluster causes these galaxy halo masses to be lower than 
that of equivalent luminosity field galaxies at comparable 
redshifts obtained from galaxy-galaxy lensing. The fraction 
of mass in these clumps is only 10-20 per cent of the total 
mass of the cluster within the inner 500 h^ 1 kpc of these high 
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Table 2. Properties for the primary and (where relevant) secondary mass clumps in the clusters. The characteristic parameters arc 
initially constrained by the positions, shapes and luminosities of the multiple-imaged objects and are then iteratively varied to match 
the weak shear field as well to obtain the optimal values (in the \ 2 sense) in a likelihood scheme. 
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Figure 1. The retrieved best-fit values for the 2 parameters that were marginalised over in the likelihood analysis, the central velocity 
dispersion and aperture radius for an L* galaxy in each of the clusters. Included is the previously published data point for AC114. The 
error bars are 3cr for both o-q* and rt* . It is these error bars that translate into a factor of 2 error in the total aperture mass. 



central density clusters. The remaining 80-90 per cent of the 
cluster mass is consistent with being smoothly distributed 
(in clumps with mass M < 10 10 Mq , the precise compo- 
sition of this component depends on the hitherto unknown 
nature of dark matter. Note that the upper and lower limits 



on the mass spectrum vary from cluster to cluster due to 
the difference in the luminosity functions of cluster galax- 
ies. These mass functions and other detailed properties of 
the substructure can now be directly compared to the sub- 
halo mass functions of dark matter haloes taken from the 
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Millennium simulation, the results of which are presented in 
Section 4. 



3.1.1 Uncertainties, sources of error and robustness of the 
lens models 

The following tests were performed for each cluster, (i) 
choosing random locations (instead of bright, early-type 
cluster member locations) for the perturbers; (ii) scrambling 
the shapes of background galaxies; (iii) choosing to asso- 
ciate the perturbers with the 40 faintest (as opposed to the 
40 brightest) galaxies; (iv) randomly selecting known cluster 
galaxies as perturbers; (v) selecting late-type galaxies. None 
of the above cases (i)-(v) yields a convergent likelihood map, 
in fact all that is seen in the resultant 2-dimensional likeli- 
hood surfaces is noise. We averaged the shear field in 4 ra- 
dial bins around each cluster galaxy and stacked for each 
cluster to estimate the aperture mass within those bins. 
We found that the mean value of the shear in the inner- 
most radial bin was ~ 8 ± 2 % falling off to ~ 2% in the 
last radial bin indicating the presence of an aperture mass 
ranging from 10 11 — 10 12 Mq. Choosing random locations in 
the cluster (i.e. not the locations of early-type cluster mem- 
bers) we do not detect a shear profile in the stacked, direct 
averaging. Alternatively, we also used the locations of late- 
type galaxies (morphological classification from HST data 
is available for all the clusters studied here), to perform the 
direct averaging and did not detect any shear. 

The robustness of our results has been extensively 
tested, however there are a couple of caveats that we ought 
to mention. As outlined above, in this galaxy-galaxy lensing 
technique we are only sensitive to a restricted mass range in 
terms of secure detection of substructure. This is due to the 
fact that we are quantifying a differential signal above the 
average tangential shear induced by a cluster, and we are 
inherently limited on by the average number of distorted 
background galaxies that lie within the aperture scale radii 
of cluster galaxies. This trade-off between requiring a suffi- 
cient number of lensed background galaxies in the vicinity of 
the subhaloes and the optimum locations for the subhaloes 
leads us to choose the brightest 40 early-type cluster galax- 
ies for each lens. With deeper, wider and more numerous 
images of clusters, expected in the future with a wide-field 
imager in space such as the SNAP mission (Aldering et al. 
2003), this technique can be pushed much further to probe 
down to lower masses in the mass function. It is possible that 
the bulk of the mass in subhaloes is in lower mass clumps 
(which in this analysis is essentially accounted for as part of 
the smooth component) and are in fact anti-correlated with 
positions of early-type galaxies. 

Our results still hold true since we are filtering out only 
the most massive clumps via this technique. Note that one 
of the null tests performed above, associating galaxy haloes 
with random positions in the cluster (and not with the lo- 
cations of bright, early-type galaxies) resulted in pure noise 
with an average value of the measured tangential shear rang- 
ing from —0.10 to 0.05. Even if we suppose that the bulk 
of the dark matter is associated with say, dwarf/very low 
surface brightness galaxies in clusters, then the spatial dis- 
tribution of these galaxies is required to be fine-tuned such 
that these effects do not show up in the shear field in the in- 
ner regions, implying that if at all they are likely to be more 



significant repositories of mass perhaps in the outskirts of 
clusters. 

Guided by the current theoretical understanding of the 
assembly of clusters, dwarf galaxies are unlikely to survive 
in the high density core regions of galaxy clusters studied 
here. Studies such as the one presented here when applied on 
larger scales to distances of a few Mpc from the cluster centre 
will enable an understanding of the role of environment in 
mass stripping of these haloes. In a recent study, Limousin 
et al. (2006) track the mass inferred for a fiducial cluster 
galaxy as a function of cluster-centric distance out to twice 
the virial radius. Using ground-based data they do not detect 
any significant variation in mass, however, using mosaiced 
HST images Natarajan et al.(2006), find that a typical L* 
galaxy-halo has higher aperture mass in the outskirts of the 
cluster. 

The principal sources of error in the above analysis are 
(i) shot noise - we are inherently limited by the finite number 
of sources sampled within a few tidal radii of each cluster 
galaxy; (ii) the spread in the intrinsic ellipticity distribution 
of the source population; (iii) observational errors arising 
from uncertainties in the measurement of ellipticities from 
the images for the faintest objects and (iv) contamination by 
foreground galaxies mistaken as background. As mentioned 
in Section 2.2.2, the partitioning of mass into subhaloes and 
the smooth component as done here is largely independent 
of the N(z) of background galaxies. 

In terms of the total contribution to the error budget, 
based on simulations we find that the shot noise is the most 
significant source of error, amounting to ~ 50 per cent; fol- 
lowed by the width of the intrinsic ellipticity distribution 
which contributes ~ 20 per cent, and the other three sources 
together contribute ~ 30 per cent. This elucidates the future 
strategy for such analyses - going significantly deeper and 
wider in terms of the field of view is likely to provide con- 
siderable gains. Mosaic-ed ACS images are the ideal data 
sets for this galaxy-galaxy lensing analyses, and such work 
is currently in progress. 



4 COMPARISON WITH N-BODY 
SIMULATIONS 

In this section, we present a comparison between the prop- 
erties of the subhaloes determined using the lensing analysis 
detailed above and results from N-body simulations. For this 
study, we make use of the Millennium Simulation, recently 
carried out by the Virgo Consortium and described in detail 
in Springel et al. (2005). The simulation follows N = 2160 3 
particles within a comoving box of size 500/i -1 Mpc on a 
side. With a particle mass of 8.6 x 10 8 fc _1 Mg and a spa- 
tial resolution of 5/i _1 kpc, the Millennium Simulation pro- 
vides an unprecedented combination of high resolution and 
large volume that allows the formation of rare objects - like 
massive lensing clusters - to be followed in a representa- 
tive fashion. It therefore represents an ideal simulation to 
compare with observations as a relatively large sample of 
massive, lensing clusters can be extracted over the redshift 
range probed by the observations. For the 64 time slices 
produced by the simulation, embedded substructures within 
dark matter haloes were identified with the algorithm SUB- 
FIND (Springel et al. 2001). In brief, the algorithm decom- 
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poses a given particle group (previously identified with a 
standard friends-of-friends algorithm) into a set of disjoint 
substructures, each of which is identified as a locally over- 
dense region in the density field of the background halo. Af- 
ter the regions containing substructure candidates have been 
identified, an unbinding procedure is employed to iteratively 
reject all particles with positive total energy. All substruc- 
tures that survive the unbinding procedure, and still have 
at least 20 self-bound particles, are considered to be genuine 
substructures. We refer to the original paper for more details 
on the algorithm. We note however that the identification 
of subhaloes within haloes represents a difficult technical 
problem. A variety of different algorithms have been devel- 
oped in order to accomplish this task (e.g. Gottlober et al. 
1998; Klypin et al. 1999, Springel et al. 2001; Weller et al. 
2005), but unfortunately little work has been done so far to 
compare the properties of subhaloes identified with differ- 
ent methods, so that the systematic effects and differences 
inherent in these methods are largely unknown. 

In order to carry out our comparison with lensing re- 
sults, we have selected all haloes with M200 > 8 x 10 14 M© 
from the simulation box at the snapshots corresponding to 
the redshifts of the lensing clusters (see Table 1). We find 
17, 15, and 12 very massive haloes at the redshifts of the 
clusters A 2218, A 2390, and CI 2244-02 respectively. Mas- 
sive cluster haloes are on the tail of the mass function in a 
ACDM Universe and their number decreases rapidly with in- 
creasing redshift. For CI 0024+16 and CI 0054-27, we have 
therefore combined two adjacent snapshots (in both cases 
the redshift of the observed clusters lies in between those 
of the used snapshots), which led to 12 and 6 candidates, 
respectively. 

For each cluster halo, we have analysed the distribution 
of the subhalo properties by projecting along the x, y, and 
z axes and keeping subhaloes 2Mpc from the cluster cen- 
tre that are projected to within 1 Mpc, identified with the 
position of the most bound particle. The distributions we 
discuss in the following are then the average over the three 
projections of all the haloes identified for each redshift. 

4.1 The mass function of substructures 

As mentioned in Sec. 1, a direct mapping of the substruc- 
ture mass function is ultimately of interest since it provides 
an important test of the ACDM paradigm. In addition, it is 
intimately connected to the galaxy formation process: com- 
parison of the dark halo mass function with the observed 
luminosity function of galaxies can provide valuable insights 
into the physical processes driving galaxy formation and 
evolution. In this work, we restrict ourselves to a compari- 
son of the properties of the dark matter component. A yet 
more detailed comparison for A 2218 in terms of observable 
galaxy properties by means of simulations and semi-analytic 
techniques is presented elsewhere (De Lucia, Natarajan & 
Springel 2007). 

In Fig. 2 we compare the substructure mass function 
retrieved from the galaxy-galaxy lensing analysis (hashed 
histograms) with the results from the numerical simulation. 
The black solid line in each panel represents the average over 
the three orthogonal projections of all the massive haloes 
identified for each redshift. All subhaloes along the line of 
sight that are projected to within 1 Mpc are included in the 



inventory. The grey filled histograms show for each value 
of the subhalo mass, the minimum and maximum number 
of substructures found in the simulated clusters. The lower 
mass cut-off for the mass function determined from lensing 
is set by observational limitations, primarily the accuracy 
with which anisotropies in the shear field can be statistically 
detected. We note that the mass of the subhaloes used to 
build the grey histogram in Fig. 2 is denned in terms of 
the total number of particles they contain, and has been 
multiplied by a correction factor of 2, as discussed below. As 
shown in Fig. 1, the error bars in ao* and rt* retrieved from 
the lensing observations translate into a factor of 2 error in 
the aperture mass of substructures. 

When taken at face value this correction factor of 2 
would signal a discrepancy between the simulation and lens- 
ing results. However, we argue that the difference is domi- 
nated by systematic effects in the measurement of substruc- 
ture masses in the simulation, and possibly also in the lens- 
ing analysis. Note that the algorithm SUBFIND identifies 
substructures as overdense regions bounded by an isoden- 
sity surface where the density equals the local value of the 
background host halo. Because of this criterion for identify- 
ing substructure, it is possible that our subhalo masses are 
biased low, because only the 'tip of the iceberg' is seen. This 
effect should increase the closer a subhalo is located to the 
centre. We will now show that this suspicion is indeed true. 

For our comparison with lensing data, the most im- 
portant question is whether the self-bound mass reported 
by SUBFIND for a substructure actually accounts for the 
full mass enhancement present at the location of the sub- 
structure in the N-body simulation. To determine this mass 
enhancement, one can measure the enclosed mass within a 
little sphere around the location of the substructure and sub- 
tract from it the mass in the same sphere after spherically 
averaging the whole mass distribution of the halo around 
the halo centre. The latter gives an estimate of the back- 
ground density in the volume occupied by the substructure. 
If the mass distribution is a superposition of a spherically 
symmetric background halo and an embedded substructure, 
then this procedure will (almost exactly) recover the true 
mass of the substructure independent of the density profile 
of the background halo and that of the substructure, pro- 
vided that the spherical aperture fully encloses the substruc- 
ture. The mass estimate will be biased slightly low because 
the spherically averaged mass distribution includes the sub- 
structure mass, but this should be negligible if the volume 
of the substructure is small against that of the background 
halo. 

This leaves the question of what to choose for the ra- 
dius of the sphere that is needed in the measurement. If 
there was only a single substructure in a spherically sym- 
metric halo, the result would be invariant once the sphere is 
large enough to fully enclose the substructure. In real-world 
haloes we will suffer from growing errors in the measure- 
ment when the aperture becomes too large, both from other 
nearby substructures and from the fact that the background 
halo is not exactly spherical in shape. As a simple compro- 
mise we have adopted 3 times the half-mass radius measured 
for the substructure (the half-mass radius is the radius that 
enclosed half the bound particles) to carry out our mea- 
surements of the mass enhancements. Note that for a NFW 
halo, the virial radius would be 2.8 times the half-mass ra- 
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Figure 2. Comparison between substructure mass function retrieved from the galaxy-galaxy lensing analysis (red shaded histograms) 
and results from haloes selected from the Millennium Simulation. The black solid line in each panel represents the average subhalo mass 
function of haloes selected at the redshift of the observed lensing cluster (see text for details). The grey shaded region represents, for 
each value of the subhalo mass, the min-max number of substructures found in the simulated haloes. 



dius for a concentration of 10, while it would be 3.2 times 
the half-mass radius for a concentration of 15. 

Adopting the above measurement procedure, we have 
determined the mass enhancements m en h around the posi- 
tions of all substructures identified by SUBFIND in haloes 
with masses above 10 14 /i _1 Mq, and with a minimum par- 
ticle number of 1000 (in total ~ 31000 objects). 

In Fig. |3] we show the ratio between the measured mass- 
enhancement and the corresponding SUBFIND self-bound 
mass as a function of radial position of the substructure 
within the parent halo. There is clearly a mass-bias with 
a strong radial dependence. The m su b masses reported by 
SUBFIND for substructures in the inner parts of clusters are 
systematically biased low; the real mass enhancements rele- 
vant for lensing can be a factor 2-3 higher in the innermost 
parts of the cluster. For the bulk of the substructures at 
larger radii, there is only a small difference of order ~ 20% 
or so, a value that could also be weakly affected by our 
choice of aperture radius. However, the radial trend of the 
mass ratios is robust. While most of the overall substruc- 
ture mass is found in the outer parts of the cluster, we note 
that the luminous cluster galaxies are much more centrally 
concentrated (Springel et al. 2001), so that the bias in the 
inner parts enters with a large weight in our comparison of 
lensing data with the N-body simulations. This provides the 
quantitative basis for our fiducial correction factor of 2 that 
we applied in the comparison shown in Fig. 2. We think 



a factor of this order is plausible given the systematics of 
current techniques, but it is also clear that a more precise 
determination will require an improved understanding of the 
systematics of both the substructure finder and the lensing 
analysis. 

We note that we have also checked whether there is any 
mass dependence of this 'bias' in the SUBFIND masses, but 
we have found no evidence for this. This is consistent with 
the fact that there is a good uniformity of the measured 
slopes for the subhalo mass function from different studies 
(De Lucia et al. 2004; Diemand et al. 2004), which suggests 
that the systematics of the substructure finders can be ac- 
counted for through simple scale factors, at least at a sta- 
tistical level. Also, we think that it is likely that other algo- 
rithms for identifying substructures suffer from this problem 
in similar ways. 

In addition to the above, is is likely that there are some 
systematic effects in the mass estimate based on the lens- 
ing analysis as well. Here a systematic overestimate appears 
more likely, arising for example if the most massive sub- 
structures correspond to haloes that have fallen in most re- 
cently, as indicated by recent numerical studies (De Lucia 
et al. 2004; Gao et al. 2004). As for the mass reconstruction 
technique employed here, biases can also be introduced due 
to the possible existence of substructure along the line of 
sight, not in the lens plane but behind it (Natarajan 2006). 
Fig. 2 also shows that there is a large scatter among simu- 
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Figure 3. Comparison of the actual mass enhancements found 
around the locations of substructures (m en h) with the self-bound 
mass m sub reported by SUBFIND, as a function of radial distance 
R s in their parent halo of virial radius R v j r . All substructures 
with more than 1000 particles in haloes above 10 14 /i _1 Mq in 
the Millennium were analysed. A random 25%-subsample of the 
points are included in the scatter diagram. The solid-line with 
symbols gives the mean in each radial bin for the whole sample, 
while the dashed line is the median in each bin. 



lated clusters and since the number of substructures in the 
observable mass range is quite small, we also expect large 
system-to-system variations between different observed sys- 
tems. Besides, as shown in Table 1, the subhalo masses are 
uncertain to within a factor of 2 from the galaxy-galaxy 
lensing analysis, based on the observational limitations from 
current data. 

We note here in our preliminary comparison of the lens- 
ing derived mass function with that extracted from the sim- 
ulation of a single, massive cluster (Natarajan & Springel 
2004), we had found excellent agreement without the need 
for a mass correction factor. However, it turns out that this 
result was influenced by not correctly taking into account the 
spatial geometry of the HST field. When this effect is appro- 
priately included, our present results are consistent with this 
earlier work, and a mass correction factor of 2 is required to 
obtain quantitative consistency. Note that this is within the 
limits of the systematic uncertainties both of the substruc- 
ture mass detection in the simulations as well as those from 
the lensing analysis. We appear to be missing subhaloes at 
the high mass end (logM > 12.5), in the lensing determi- 
nation. This is likely due to the fact that such high mass 
haloes probably host more than one luminous component. 
We plan to investigate this issue further in our detailed study 
of A2218 in future work. 

Overall, it is therefore remarkable that we find very 
good qualitative agreement in the mass range sampled by 
the lensing analysis, and very good quantitative agreement 
as well when the systematic biases are approximately cor- 



rected for within the errors. It is interesting to note here 
that the only obvious outlier is CI 0024+16. This is the only 
system that clearly has a bimodal structure with two almost 
equal mass clumps separated by a small redshift offset. Its 
bimodal structure has been interpreted as a sign of a re- 
cent merger (Czoske et al. 2002). Clearly, if this is indeed 
the case, the comparison with relaxed massive haloes from 
the Millennium Simulation is not appropriate. We plan to 
investigate this question in more detail in future work by 
comparing the properties of the substructures in CI 0024+16 
with clusters selected from the Millennium Simulation that 
are recent mergers of equivalent mass. 



4.2 The distribution of velocity dispersions 

In Fig. 4, we compare the distribution of velocity disper- 
sions retrieved by the lensing analysis and the average dis- 
tributions obtained for the simulated clusters selected from 
the Millennium Simulation. The velocity dispersions of sub- 
structures in simulated clusters are estimated using the ve- 
locity information for all the particles attached to the dark 
matter substructure and have been scaled by a factor of 
v2 (see previous section). Again we find generally a good 
agreement over the range sampled by the lensing analysis, 
except for CI 0024+16, which yet again appears to be an 
outlier (see discussion in previous section). Note here that 
the velocity dispersion retrieved from the lensing analysis 
is one of the two parameters that characterises each sub- 
halo. This velocity dispersion is the normalization of the 
Faber- Jackson relation, therefore, relates the mass (for the 
PIEMD model) to that of early-type galaxies. For the simu- 
lated haloes the measured velocity dispersion is that of the 
dark matter alone. We note that Springel et al. (2001) used 
a factor equal to 0.9 to convert the ID dispersions of sub- 
haloes into stellar velocity dispersions, which was adopted 
to match the zero-point of the Faber-Jackson relation. Due 
to the dissipation involved in star formation, it is plausible 
that the stellar velocity dispersion is smaller than that of 
the dark matter. On the other hand, we argued above that 
our mass estimates might be biased low. We explore spatial 
and velocity biases in more detail using the lensing results 
from A 2218 in a forthcoming paper. 



4.3 The distribution function of tidal radii 

We use the aperture radii derived from the lensing analy- 
sis as a proxy for the tidal radii for the subhaloes. These 
aperture radii are more compact than the tidally truncated 
extents of the subhaloes. This radius cannot be identified 
directly with the tidal radius defined in the usual sense, but 
is proportional to it. In Fig. 5, we finally show the com- 
parison between the distribution of aperture radii (rt) re- 
trieved by the lensing analysis with the average distribu- 
tions of the tidal radii of substructures in simulated clusters 
selected from the Millennium Simulation. The solid black 
line in the figure and the grey filled region are obtained by 
measuring, for each substructure, the radius rhaif that di- 
vides the halo mass in equal inner and outer parts. We note 
that, by using this measure of the radius, we find a very 
nice agreement with the distributions retrieved by the lens- 
ing analysis (with the usual exception of CI 0024+16) over 
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Figure 4. As in Fig. 2 but for the distribution of velocity dispersions. 



the full range sampled by the observations. ri la if, however, 
represents an underestimate of the 'tidal radius' of a dark 
matter substructure. The dashed black line in Fig. 5 gives 
the distribution obtained by using the radius corresponding 
to the maximum circular velocity (r max ), which is offset with 
respect to rhaif by a factor that is slightly less than 2. We 
note that, as explained in the previous sections, the mass re- 
trieved by the lensing analysis is essentially the mass within 
a fixed aperture, that is then identified as being proportional 
to the tidal radius of the associated dark matter substruc- 
ture. Given the uncertainty in the measurement technique, 
the good agreement shown in Fig. 5 is therefore notewor- 
thy. We recall that, in order to obtain a good agreement 
with the observed subhalo mass function, we have corrected 
their mass by a factor of 2. We are therefore in a situa- 
tion where the substructures detected in the simulations are 
slightly bigger but less massive than those retrieved by our 
lensing technique. Such a situation would be expected if the 
density profiles assumed for the substructures are systemati- 
cally denser in the inner parts than the simulated subhaloes. 
And, in fact, recent numerical studies suggest that dark mat- 
ter subhaloes have softer profiles than NFW (Hayashi et al. 
2004; Stoehr et al. 2003). Unfortunately our lensing tech- 
nique is not sensitive to the choice of mass profile, rather 
it represents a tool to determine the total mass enclosed 
within an aperture. Higher resolution data will in the future 
allow the slopes of the mass profiles in substructure to be 
constrained. Adiabatic contraction of the dark matter den- 
sity profile in response to the baryonic component in the 
inner regions of the cluster halo (on scales of r < 10 — 50 



kpc) is also expected to be relevant, however the importance 
and significance remains to be assessed from numerical sim- 
ulations. There are two conflicting claims on how significant 
this effect at the present time (Zappacosta et al. 2006; Rozo, 
Nagai, Keeton & Kravtsov 2006). 



5 CONCLUSIONS AND DISCUSSION 

In this paper, we present (i) high resolution mass models for 
lensing clusters, (ii) the inferred mass function of subhaloes 
inside these clusters and (iii) a detailed comparison of the 
subhalo mass function, the velocity dispersion and aperture 
radii function with an ensemble of cluster-sized haloes se- 
lected from the Millennium Simulation. Detailed results of 
the application of our galaxy-galaxy lensing analysis tech- 
niques to five HST cluster lenses are used to construct high 
resolution mass models of the inner regions. In order to do so 
we have utilised both strong and weak lensing observations 
for these massive clusters. The goal has been to quantify sub- 
structure in the cluster assuming that the subhaloes follow 
the distribution of bright, early-type cluster galaxies. Simi- 
lar attempts have been made in the lower density field envi- 
ronment yielding typical galaxy masses and central velocity 
dispersions. The mass distribution for a typical galaxy halo 
inferred from field studies is extended with no discernible 
cut-off in most analyses with the exception of one study by 
Hoekstra et al.(2004) where they report r t = 260t 7 3 4 kpc 
using data from the CNOC-2 fields. By contrast, in the clus- 
ter environments probed in this work we detect an edge to 
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Figure 5. As in Fig. 2 but for the distribution of tidal radii. The solid black line with the grey regions are obtained by measuring for 
each substructure, the radius corresponding to its half mass. The dashed black line is obtained by taking the radius corresponding to 
the maximum circular velocity. 



the mass distribution in cluster galaxies. We have performed 
various stringent checks to ascertain that this is not an arti- 
fact of the choice of mass model. Rather this result provides 
evidence for tidal stripping by the global cluster potential 
well. 

Aside from the detailed lens models, we also present 
the mass spectrum (albeit within a limited mass range with 



subhalo masses ranging from 10 



10 1 



Mq) of substruc- 



ture in the inner regions of these clusters. The survival and 
evolution of substructure offers a stringent test of structure 
formation models within the CDM paradigm. Subhaloes of 
the scale detected in all these clusters indicate a significant 
probability of galaxy-galaxy collisions over a Hubble time 
within a rich cluster. However, since the internal velocity 
dispersions of these clumps associated with early-type clus- 
ter galaxies (~ 150-250 kms ) are much smaller than their 
orbital velocities, these interactions are unlikely to lead to 
mergers, suggesting that the encounters of the kind simu- 
lated in the galaxy harassment picture by Moore et al. (1996) 
are the most frequent and likely. High resolution cosmolog- 
ical N-body simulations of cluster formation and evolution 
(De Lucia et al. 2004; Ghigna et al. 1998; Moore et al. 1996), 
find that the dominant interactions are between the global 
cluster tidal field and individual galaxies after z — 2. The 
cluster tidal field tidally strips galaxy haloes in the inner 0.5 
Mpc efficiently and the radial extent of the surviving haloes 
is a strong function of their distance from the cluster centre. 



Much of this modification is found to occur between z = 
and z = 0.5. 

We have performed a detailed comparison of the sub- 
halo mass function, velocity dispersion function and the dis- 
tribution of aperture radii retrieved from lensing with an 
ensemble of massive clusters from the Millennium Simula- 
tion, which provides an ideal data-set for such an analysis. 
A series of systematic effects, due to the algorithm used to 
identify substructures in the simulation, needs to be taken 
into account when performing such a comparison. In addi- 
tion, it should be kept in mind that the level of uncertainty of 
the observational results is still relatively high (for instance 
the mass derived from lensing is only accurate to within a 
factor of two and the technique is not sensitive to the choice 
of the substructure mass profile). Overall, we find consis- 
tency between the distribution of substructure properties 
retrieved using the lensing analysis and those obtained from 
simulations, although our detailed comparison seems to sug- 
gest that simulated substructures are slightly bigger but less 
massive than sub-clumps detected by means of lensing tech- 
niques. This might be due to systematic differences between 
the density profiles of simulated substructures and those as- 
sumed in the lensing model. Unfortunately, the technique is 
not sensitive to this choice but higher resolution data will 
allow in the future the slopes of the mass profiles in sub- 
structures to be constrained. 

Despite the uncertainties mentioned above, the general 
agreement between simulations and results determined inde- 
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pendently from lensing is remarkable. Our work represents a 
powerful test of the ACDM model, which at present appears 
to be consistent with the amount of observed substructure 
in massive, lensing clusters, up to redshifts of ~ 0.6, given 
the uncertainties. It will be very interesting to tighten the 
constraints with future lensing data. 
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